Library import
library(SpatialExperiment)
library(data.table)
library(scater) # it imports scuttle
library(cowplot)
library(ggplot2)
theme_set(theme_bw())
library(matrixStats)
library(dplyr)
library(tidyr)
library(tibble)
library(arrow)
library(scales)
library(sf)
library(tmap)
library(knitr)
library(grid)
library(ggpointdensity)
library(ggpubr)
library(spdep)
library(ggExtra)
library(terra)
library(DT)
library(e1071)
library(robustbase)
library(InSituType)
library(UpSetR)
library(sfheaders)
spe <- readRDS("/data01/Shared/Benedetta/spatial_transcriptomics/CosMx/Datasets/CosMx_Breast_DBKERO/no_removal_of_outliers/no_removal_spe.rds")
sample_name <- "CosMx_DBKero_BC"
celltype_palette <- c(
"0_counts" ="black",
"APC TAMs" = "royalblue1",
"M2 TAMs" = "blue",
"cDCs" = "dodgerblue4",
"mDCs" = "deepskyblue",
"Low-quality TAMs" = "cyan",
"Mast cells" = "lightblue1",
"CTLs" = "darkseagreen1",
"NK cells" = "darkolivegreen1",
"Tregs" = "limegreen",
"Mural cells" = "goldenrod4",
"Myoepithelial cells" = "goldenrod1",
"Blood ECs" = "lightgoldenrod1",
"CAFs" = "lightgoldenrod4",
"MG Luminal cells" = "darkorange",
"Plasma cells" = "deeppink1",
"Mix BC cells TAMs" = "purple",
"ER+HER2+BC cells" = "orangered1",
"Cycling ER+HER2+BC cells" = "red4",
"Basal-like BC cells" = "darkorange3"
)
spe$polygons$control_sum <- spe$control_sum
spe$polygons$ctrl_total_ratio <- spe$ctrl_total_ratio
spe$polygons$Area_um <- spe$Area_um
spe$polygons$log2AspectRatio <- spe$log2AspectRatio
spe$polygons$Mean.DAPI <- spe$Mean.DAPI
spe$target_area_ratio <- spe$target_sum/spe$Area_um
spe$polygons$target_area_ratio<- spe$target_area_ratio
To view the spatial organization of both the sample as it is and FOVs
metadata <- data.frame(colData(spe))
# image theme
global_map_theme <- function(back.color="black",back.border=NA,title.col="white") {
theme(aspect.ratio = 1,
panel.border = element_blank(),
legend.key = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.ticks = element_blank(),
axis.text.y = element_blank(),
axis.text.x = element_blank(),
panel.grid = element_blank(),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_rect(fill = "transparent",colour = NA),
plot.title = element_text(color = title.col,hjust = 0.5,face = "bold"),
plot.background = element_rect(fill = "transparent",colour = back.border))
}
ggplot() +
geom_point(data = metadata, mapping = aes(x = CenterX_global_px, y = CenterY_global_px), colour = "darkmagenta", fill="darkmagenta", size = 0.05, alpha = 0.2) +
annotate("rect", xmin = metadata(spe)$fov_positions$x_global_px,
xmax = metadata(spe)$fov_positions$x_global_px + metadata(spe)$fov_dim["xdim"],
ymin = metadata(spe)$fov_positions$y_global_px,
ymax = metadata(spe)$fov_positions$y_global_px + metadata(spe)$fov_dim["ydim"],
alpha = .2, color = "black",linewidth = 0.2) +
geom_text(aes(x = metadata(spe)$fov_positions$x_global_px+metadata(spe)$fov_dim["xdim"]/2,
y = metadata(spe)$fov_positions$y_global_px+metadata(spe)$fov_dim["ydim"]/2,
label = metadata(spe)$fov_positions$fov), color="black") +
ggtitle(sample_name) +
global_map_theme(back.color = "white", back.border = "white", title.col = "black")
tm_shape(spe$polygons) +
tm_borders(lwd = 0.1, col = "grey50") +
tm_layout(legend.outside = TRUE,
main.title.position = c("center", "top"),
main.title = sample_name,
main.title.size = 0.5,
inner.margins = c(0, 0, 0, 0),
outer.margins = c(0, 0, 0, 0))
## Warning: The projection of the shape object spe$polygons is not known, while it
## seems to be projected.
## Warning: Current projection of shape spe$polygons unknown and cannot be
## determined.
To spot areas with anomalous feature values and spatial patterns on a global scale
plotColData(spe, "CenterY_global_um", "CenterX_global_um", colour_by = "total", point_size = 0.1) + ggtitle("Total counts per cell") + theme(aspect.ratio = 1)
plotColData(spe, "CenterY_global_um", "CenterX_global_um", colour_by = "detected", point_size = 0.1) + ggtitle("Detected features per cell") + theme(aspect.ratio = 1)
plotColData(spe, "CenterY_global_um", "CenterX_global_um", colour_by = "Area_um", point_size = 0.1) + ggtitle("Area in μm per cell") + theme(aspect.ratio = 1)
plotColData(spe, "CenterY_global_um", "CenterX_global_um", colour_by = "log2AspectRatio", point_size = 0.1) + ggtitle("Logged width/height ratio per cell") + theme(aspect.ratio = 1)
plotColData(spe, "CenterY_global_um", "CenterX_global_um", order_by = "control_sum", colour_by = "control_sum", point_size = 0.1) +
scale_color_gradient(low = "white", high = "red", name = "control_sum") + ggtitle("Negative control probe counts per cell") +
theme(aspect.ratio = 1,
panel.background = element_rect(fill = "black",color = NA),
plot.background = element_rect(fill = "black",color = NA),
axis.title.x = element_text(color = "white"),
axis.title.y = element_text(color = "white"),
axis.ticks = element_line(color = "white"),
axis.text.y = element_text(color = "white"),
axis.text.x = element_text(color = "white"),
axis.line = element_line(color = "white"),
legend.title = element_text(color = "white"),
legend.text = element_text(color = "white"),
plot.title = element_text(color = "white"))
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
plotColData(spe, "CenterY_global_um", "CenterX_global_um", colour_by = "ctrl_total_ratio", point_size = 0.1) + ggtitle("Control counts/total counts ratio per cell") + theme(aspect.ratio = 1)
To check feature distribution before highlighting where extreme values are located in the global sample map
variables <- c("total", "detected", "Area_um", "log2AspectRatio", "ctrl_total_ratio", "control_sum", "Mean.DAPI")
for (var in 1:length(variables)){
n <- grep(paste0("^", variables[var], "$"), colnames(spe@colData))
histo_var <- names(colData(spe))[n]
p <- ggplot(data = data.frame(colData(spe))) +
geom_histogram(aes(x = .data[[histo_var]]), fill = "#69b3a2") +
ggtitle(paste0(histo_var))
print(p)
}
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
tm_shape(spe$polygons) +
tm_fill(col = "flagged_0total") +
tm_layout(legend.outside = TRUE,
main.title.position = c("center", "top"),
main.title = paste0("Cells having 0 total counts", "\nFlagged cells = ",
table(colData(spe)$flagged_0total)["TRUE"], " out of ", dim(spe)[2], ", ",
round(table(colData(spe)$flagged_0total)["TRUE"]/dim(spe)[2]*100, 2), "% of total cells"),
main.title.fontface = 2,
main.title.size = 1,
inner.margins = c(0, 0, 0, 0),
outer.margins = c(0, 0, 0, 0))
## Warning: The projection of the shape object spe$polygons is not known, while it
## seems to be projected.
## Warning: Current projection of shape spe$polygons unknown and cannot be
## determined.
tm_shape(spe$polygons) +
tm_fill(col= "outlier_ctrl_tot") +
tm_layout(legend.outside = TRUE,
main.title.position = c("center", "top"),
main.title = paste0("Outlier cells for control / total counts ratio > 0.1, \nflagged cells = ",
table(colData(spe)$outlier_ctrl_tot)["TRUE"], " out of ", dim(spe)[2], ", ",
round(table(colData(spe)$outlier_ctrl_tot)["TRUE"]/dim(spe)[2]*100, 2), "% of total cells"),
main.title.fontface = 2,
main.title.size = 1,
inner.margins = c(0, 0, 0, 0),
outer.margins = c(0, 0, 0, 0))
## Warning: The projection of the shape object spe$polygons is not known, while it
## seems to be projected.
## Warning: Current projection of shape spe$polygons unknown and cannot be
## determined.
out <- adjbox(spe$Area_um, plot = FALSE)
metadata(spe)$area_fence <- out$fence
tm_shape(spe$polygons) +
tm_fill(col= "umarea_outlier") +
tm_layout(legend.outside = TRUE,
main.title.position = c("center", "top"),
main.title = paste0("Outlier cells for um area",
"\n fence = " , round(metadata(spe)$area_fence[1,1], 2), ", ", round(metadata(spe)$area_fence[2,1], 2),
"\nFlagged cells = ",
table(spe$polygons$umarea_outlier)["red"], " out of ", dim(spe)[2], ", ",
round(table(spe$polygons$umarea_outlier)["red"]/dim(spe)[2]*100, 2), "% of total cells"),
main.title.fontface = 2,
main.title.size = 1,
inner.margins = c(0, 0, 0, 0),
outer.margins = c(0, 0, 0, 0))
## Warning: The projection of the shape object spe$polygons is not known, while it
## seems to be projected.
## Warning: Current projection of shape spe$polygons unknown and cannot be
## determined.
out <- adjbox(spe$Mean.DAPI, plot = FALSE)
metadata(spe)$dapi_fence <- out$fence
tm_shape(spe$polygons) +
tm_fill(col= "dapi_outlier") +
tm_layout(legend.outside = TRUE,
main.title.position = c("center", "top"),
main.title = paste0("Outlier cells for dapi",
"\n fence = " , round(metadata(spe)$dapi_fence[1,1], 2), ", ", round(metadata(spe)$dapi_fence[2,1], 2),
"\nFlagged cells = ",
table(spe$polygons$dapi_outlier)["red"], " out of ", dim(spe)[2], ", ",
round(table(spe$polygons$dapi_outlier)["red"]/dim(spe)[2]*100, 2), "% of total cells"),
main.title.fontface = 2,
main.title.size = 1,
inner.margins = c(0, 0, 0, 0),
outer.margins = c(0, 0, 0, 0))
## Warning: The projection of the shape object spe$polygons is not known, while it
## seems to be projected.
## Warning: Current projection of shape spe$polygons unknown and cannot be
## determined.
tm_shape(spe$polygons) +
tm_fill(col= "fscore_outlier") +
tm_layout(legend.outside = TRUE,
main.title.position = c("center", "top"),
main.title = paste0("Outlier cells for flag score",
"\n threshold = " , 0.6,
"\nFlagged cells = ",
table(spe$polygons$fscore_outlier)["red"], " out of ", dim(spe)[2], ", ",
round(table(spe$polygons$fscore_outlier)["red"]/dim(spe)[2]*100, 2), "% of total cells"),
main.title.fontface = 2,
main.title.size = 1,
inner.margins = c(0, 0, 0, 0),
outer.margins = c(0, 0, 0, 0))
## Warning: The projection of the shape object spe$polygons is not known, while it
## seems to be projected.
## Warning: Current projection of shape spe$polygons unknown and cannot be
## determined.
df <- data.frame(colData(spe))
colors <- c("Lower thr." = "purple4", "Upper thr." = "tomato")
ggplot(data = df) +
geom_histogram(aes(x = Area_um), fill = "#81ccbb") +
geom_vline(aes(xintercept = round(metadata(spe)$area_fence, 2)[1], color = "Lower thr.")) +
geom_vline(aes(xintercept = round(metadata(spe)$area_fence, 2)[2], color = "Upper thr.")) +
geom_text(aes(x = round(metadata(spe)$area_fence, 2)[1], y = -2083, label = as.character(round(metadata(spe)$area_fence, 2)[1])), color = "purple4") +
geom_text(aes(x = round(metadata(spe)$area_fence, 2)[2], y = -2083, label = as.character(round(metadata(spe)$area_fence, 2)[2])), color = "tomato") +
ggtitle("Cell area distribution with lower and upper threshold for outliers") +
labs(color = "Legend") +
scale_color_manual(values = colors)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(data = df) +
geom_histogram(aes(x = Mean.DAPI), fill = "#81ccbb") +
geom_vline(aes(xintercept = round(metadata(spe)$dapi_fence, 2)[1], color = "Lower thr.")) +
geom_vline(aes(xintercept = round(metadata(spe)$dapi_fence, 2)[2], color = "Upper thr.")) +
geom_text(aes(x = round(metadata(spe)$dapi_fence, 2)[1], y = -2083, label = as.character(round(metadata(spe)$dapi_fence, 2)[1])), color = "purple4") +
geom_text(aes(x = round(metadata(spe)$dapi_fence, 2)[2], y = -2083, label = as.character(round(metadata(spe)$dapi_fence, 2)[2])), color = "tomato") +
ggtitle("Mean DAPI signal with lower and upper threshold for outliers") +
labs(color = "Legend") +
scale_color_manual(values = colors)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#df <- data.frame(colData(spe))
ggplot(df, aes(x=InSituType_Refined, y = Area_um, fill = InSituType_Refined))+
geom_violin() +
scale_color_manual(values = colors)+
geom_hline(aes(yintercept = round(metadata(spe)$area_fence, 2)[1], color = "Lower thr.")) +
geom_hline(aes(yintercept = round(metadata(spe)$area_fence, 2)[2], color = "Upper thr.")) +
geom_text(aes(x = 1.5, y = round(metadata(spe)$area_fence, 2)[1] - 25,label = as.character(round(metadata(spe)$area_fence, 2)[1])), color = "purple4")+
geom_text(aes(x = 1.5, y = round(metadata(spe)$area_fence, 2)[2] + 25,label = as.character(round(metadata(spe)$area_fence, 2)[2])), color = "tomato") +
#geom_jitter(shape=16, position=position_jitter(0.2), alpha = 0.02) +
scale_fill_manual(values = celltype_palette)+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
labs(title = "Distribution of area in um values by IST cell types")
#df <- data.frame(colData(spe))
ggplot(df, aes(x=InSituType_Refined, y = Mean.DAPI, fill = InSituType_Refined))+
geom_violin() +
scale_color_manual(values = colors)+
geom_hline(aes(yintercept = round(metadata(spe)$dapi_fence, 2)[1], color = "Lower thr.")) +
geom_hline(aes(yintercept = round(metadata(spe)$dapi_fence, 2)[2], color = "Upper thr.")) +
geom_text(aes(x = 1.5, y = round(metadata(spe)$dapi_fence, 2)[1] - 250,label = as.character(round(metadata(spe)$dapi_fence, 2)[1])), color = "purple4")+
geom_text(aes(x = 2, y = round(metadata(spe)$dapi_fence, 2)[2] + 250,label = as.character(round(metadata(spe)$dapi_fence, 2)[2])), color = "tomato") +
#geom_jitter(shape=16, position=position_jitter(0.2), alpha = 0.02) +
scale_fill_manual(values = celltype_palette)+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
labs(title = "Distribution of DAPI signal values by IST cell types")
cell_df <- data.frame(colData(spe))
p1 <- ggplot(data = cell_df) +
geom_pointdensity(aes(x = log2AspectRatio, y = total)) +
scale_color_viridis_c() +
theme(legend.position = "none",
plot.margin = margin())
p2 <- ggdensity(data = cell_df, x = "total", color = "#ccccff",
fill = "#ccccff", alpha = 0.5) + theme_classic() +
theme(legend.position = "none", axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.line.y = element_blank(),
axis.line.x = element_blank(),
plot.margin = margin()) + coord_flip()
p3 <- ggdensity(data = cell_df, x = "log2AspectRatio", color = "#ccccff",
fill = "#ccccff", alpha = 0.5) + theme_classic() +
theme(legend.position = "none", axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.line.x = element_blank(),
axis.line.y = element_blank(),
plot.margin = margin())
plot_grid(p3, NULL, p1, p2, ncol = 2, align = "hv",
rel_widths = c(2, 1), rel_heights = c(1, 2))
Detected features ~ logged width/height ratio
p1 <- ggplot(data = cell_df) +
geom_pointdensity(aes(x = log2AspectRatio, y = detected)) +
scale_color_viridis_c() +
theme(legend.position = "none",
plot.margin = margin())
p2 <- ggdensity(data = cell_df, x = "detected", color = "#ccccff",
fill = "#ccccff", alpha = 0.5) + theme_classic() +
theme(legend.position = "none", axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.line.y = element_blank(),
axis.line.x = element_blank(),
plot.margin = margin()) + coord_flip()
p3 <- ggdensity(data = cell_df, x = "log2AspectRatio", color = "#ccccff",
fill = "#ccccff", alpha = 0.5) + theme_classic() +
theme(legend.position = "none", axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.line.x = element_blank(),
axis.line.y = element_blank(),
plot.margin = margin())
plot_grid(p3, NULL, p1, p2, ncol = 2, align = "hv",
rel_widths = c(2, 1), rel_heights = c(1, 2))
Total counts ~ area in um
p1 <- ggplot(data = cell_df) +
geom_pointdensity(aes(x = Area_um, y = total)) +
scale_color_viridis_c() +
theme(legend.position = "none",
plot.margin = margin())
p2 <- ggdensity(data = cell_df, x = "total", color = "#ccccff",
fill = "#ccccff", alpha = 0.5) + theme_classic() +
theme(legend.position = "none", axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.line.y = element_blank(),
axis.line.x = element_blank(),
plot.margin = margin()) + coord_flip()
p3 <- ggdensity(data = cell_df, x = "Area_um", color = "#ccccff",
fill = "#ccccff", alpha = 0.5) + theme_classic() +
theme(legend.position = "none", axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.line.x = element_blank(),
axis.line.y = element_blank(),
plot.margin = margin())
plot_grid(p3, NULL, p1, p2, ncol = 2, align = "hv",
rel_widths = c(2,1), rel_heights = c(1, 2))
Logged width/height ratio ~ negative control probe counts
p1 <- ggplot(data = cell_df) +
geom_pointdensity(aes(x = subsets_NegPrb_sum, y = log2AspectRatio)) +
scale_color_viridis_c() +
theme(legend.position = "none",
plot.margin = margin())
p2 <- ggdensity(data = cell_df, x = "log2AspectRatio", color = "#ccccff",
fill = "#ccccff", alpha = 0.5) + theme_classic() +
theme(legend.position = "none", axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.line.y = element_blank(),
axis.line.x = element_blank(),
plot.margin = margin()) + coord_flip()
p3 <- ggdensity(data = cell_df, x = "subsets_NegPrb_sum", color = "#ccccff",
fill = "#ccccff", alpha = 0.5) + theme_classic() +
theme(legend.position = "none", axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.line.x = element_blank(),
axis.line.y = element_blank(),
plot.margin = margin())
plot_grid(p3, NULL, p1, p2, ncol = 2, align = "hv",
rel_widths = c(2,1), rel_heights = c(1, 2))
Logged width/height ratio ~ local coordinates on X axis in px
p1 <- ggplot(data = cell_df) +
geom_pointdensity(aes(x = CenterX_local_px, y = log2AspectRatio)) +
scale_color_viridis_c() +
theme(legend.position = "none",
plot.margin = margin())
p2 <- ggdensity(data = cell_df, x = "log2AspectRatio", color = "#ccccff",
fill = "#ccccff", alpha = 0.5) + theme_classic() +
theme(legend.position = "none", axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.line.y = element_blank(),
axis.line.x = element_blank(),
plot.margin = margin()) + coord_flip()
p3 <- ggdensity(data = cell_df, x = "CenterX_local_px", color = "#ccccff",
fill = "#ccccff", alpha = 0.5) + theme_classic() +
theme(legend.position = "none", axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.line.x = element_blank(),
axis.line.y = element_blank(),
plot.margin = margin())
plot_grid(p3, NULL, p1, p2, ncol = 2, align = "hv",
rel_widths = c(2,1), rel_heights = c(1, 2))
Logged width/height ratio ~ local coordinates on X axis in px
p1 <- ggplot(data = cell_df) +
geom_pointdensity(aes(x = CenterY_local_px, y = log2AspectRatio)) +
scale_color_viridis_c() +
theme(legend.position = "none",
plot.margin = margin())
p2 <- ggdensity(data = cell_df, x = "log2AspectRatio", color = "#ccccff",
fill = "#ccccff", alpha = 0.5) + theme_classic() +
theme(legend.position = "none", axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.line.y = element_blank(),
axis.line.x = element_blank(),
plot.margin = margin()) + coord_flip()
p3 <- ggdensity(data = cell_df, x = "CenterY_local_px", color = "#ccccff",
fill = "#ccccff", alpha = 0.5) + theme_classic() +
theme(legend.position = "none", axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.line.x = element_blank(),
axis.line.y = element_blank(),
plot.margin = margin())
plot_grid(p3, NULL, p1, p2, ncol = 2, align = "hv",
rel_widths = c(2,1), rel_heights = c(1, 2))
spe$polygons$fov_id <- paste(spe$polygons$fov, spe$polygons$cellID, sep = "_")
spe$polygons$outlier_fence_color <- case_when(spe$flagged_0total == TRUE ~ "darkturquoise",
spe$outlier_ctrl_tot == TRUE ~ "magenta",
spe$Mean.DAPI > round(metadata(spe)$dapi_fence[2,1], 2) ~ "greenyellow",
spe$Mean.DAPI < round(metadata(spe)$dapi_fence[1,1], 2) ~ "purple",
spe$Area_um > round(metadata(spe)$area_fence[2,1], 2) ~ "red",
spe$Area_um < round(metadata(spe)$area_fence[1,1], 2) ~ "white",
TRUE ~ "black")
polygons_sub <- spe$polygons[spe$polygons$fov%in%1:8,]
tm_shape(polygons_sub) +
tm_fill(col= "outlier_fence_color") +
tm_borders(lwd = 0.2, col = "cyan") +
tm_layout(legend.outside = TRUE,
main.title.position = c("center", "top"),
main.title = "",
main.title.fontface = 2,
main.title.size = 1,
inner.margins = c(0, 0, 0, 0),
outer.margins = c(0, 0, 0, 0),
bg.color = "black") +
tm_add_legend(col = c("black", "darkturquoise", "magenta", "purple", "greenyellow","white", "red"),
labels = c("unflagged cells", "cells with 0 total counts", "cells with ctrl/total ratio > 0.1", "outlier for DAPI lower thr.", "outlier for DAPI higher thr.", "outlier for um area lower thr.","outlier for um area higher thr."))
## Warning: The projection of the shape object polygons_sub is not known, while it
## seems to be projected.
## Warning: Current projection of shape polygons_sub unknown and cannot be
## determined.
## Some legend labels were too wide. These labels have been resized to 0.62, 0.65, 0.63. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
tm_shape(polygons_sub) +
tm_fill(col= "IST_Refined_color") +
tm_layout(legend.outside = TRUE,
main.title.position = c("center", "top"),
main.title = "",
main.title.fontface = 2,
main.title.size = 1,
inner.margins = c(0, 0, 0, 0),
outer.margins = c(0, 0, 0, 0),
bg.color = "black") +
tm_add_legend(labels = names(celltype_palette),
col = celltype_palette)
## Warning: The projection of the shape object polygons_sub is not known, while it
## seems to be projected.
## Warning: Current projection of shape polygons_sub unknown and cannot be
## determined.
polygons1 <- subset(polygons_sub, subset = fov == 1)
tm_shape(polygons1) +
tm_fill(col= "outlier_fence_color") +
tm_borders(lwd = 0.2, col = "cyan") +
tm_layout(legend.outside = TRUE,
main.title.position = c("center", "top"),
main.title = "",
main.title.fontface = 2,
main.title.size = 1,
inner.margins = c(0, 0, 0, 0),
outer.margins = c(0, 0, 0, 0),
bg.color = "black")+
tm_add_legend(col = c("black", "darkturquoise", "magenta", "purple", "greenyellow","white", "red"),
labels = c("unflagged cells", "cells with 0 total counts", "cells with ctrl/total ratio > 0.1", "outlier for DAPI lower thr.", "outlier for DAPI higher thr.", "outlier for um area lower thr.","outlier for um area higher thr."))
## Warning: The projection of the shape object polygons1 is not known, while it
## seems to be projected.
## Warning: Current projection of shape polygons1 unknown and cannot be
## determined.
## Some legend labels were too wide. These labels have been resized to 0.62, 0.65, 0.63. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
polygons2 <- subset(polygons_sub, subset = fov == 2)
tm_shape(polygons2) +
tm_fill(col= "outlier_fence_color") +
tm_borders(lwd = 0.2, col = "cyan") +
tm_layout(legend.outside = TRUE,
main.title.position = c("center", "top"),
main.title = "",
main.title.fontface = 2,
main.title.size = 1,
inner.margins = c(0, 0, 0, 0),
outer.margins = c(0, 0, 0, 0),
bg.color = "black")+
tm_add_legend(col = c("black", "darkturquoise", "magenta", "purple", "greenyellow","white", "red"),
labels = c("unflagged cells", "cells with 0 total counts", "cells with ctrl/total ratio > 0.1", "outlier for DAPI lower thr.", "outlier for DAPI higher thr.", "outlier for um area lower thr.","outlier for um area higher thr."))
## Warning: The projection of the shape object polygons2 is not known, while it
## seems to be projected.
## Warning: Current projection of shape polygons2 unknown and cannot be
## determined.
## Some legend labels were too wide. These labels have been resized to 0.62, 0.65, 0.63. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
polygons3 <- subset(polygons_sub, subset = fov == 3)
tm_shape(polygons3) +
tm_fill(col= "outlier_fence_color") +
tm_borders(lwd = 0.2, col = "cyan") +
tm_layout(legend.outside = TRUE,
main.title.position = c("center", "top"),
main.title = "",
main.title.fontface = 2,
main.title.size = 1,
inner.margins = c(0, 0, 0, 0),
outer.margins = c(0, 0, 0, 0),
bg.color = "black")+
tm_add_legend(col = c("black", "darkturquoise", "magenta", "purple", "greenyellow","white", "red"),
labels = c("unflagged cells", "cells with 0 total counts", "cells with ctrl/total ratio > 0.1", "outlier for DAPI lower thr.", "outlier for DAPI higher thr.", "outlier for um area lower thr.","outlier for um area higher thr."))
## Warning: The projection of the shape object polygons3 is not known, while it
## seems to be projected.
## Warning: Current projection of shape polygons3 unknown and cannot be
## determined.
## Some legend labels were too wide. These labels have been resized to 0.62, 0.65, 0.63. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
polygons3 <- subset(polygons_sub, subset = fov == 3)
tm_shape(polygons3) +
tm_fill(col= "outlier_fence_color") +
tm_borders(lwd = 0.2, col = "cyan") +
tm_layout(legend.outside = TRUE,
main.title.position = c("center", "top"),
main.title = "",
main.title.fontface = 2,
main.title.size = 1,
inner.margins = c(0, 0, 0, 0),
outer.margins = c(0, 0, 0, 0),
bg.color = "black")+
tm_add_legend(col = c("black", "darkturquoise", "magenta", "purple", "greenyellow","white", "red"),
labels = c("unflagged cells", "cells with 0 total counts", "cells with ctrl/total ratio > 0.1", "outlier for DAPI lower thr.", "outlier for DAPI higher thr.", "outlier for um area lower thr.","outlier for um area higher thr."))
## Warning: The projection of the shape object polygons3 is not known, while it
## seems to be projected.
## Warning: Current projection of shape polygons3 unknown and cannot be
## determined.
## Some legend labels were too wide. These labels have been resized to 0.62, 0.65, 0.63. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
polygons4 <- subset(polygons_sub, subset = fov == 4)
tm_shape(polygons4) +
tm_fill(col= "outlier_fence_color") +
tm_borders(lwd = 0.2, col = "cyan") +
tm_layout(legend.outside = TRUE,
main.title.position = c("center", "top"),
main.title = "",
main.title.fontface = 2,
main.title.size = 1,
inner.margins = c(0, 0, 0, 0),
outer.margins = c(0, 0, 0, 0),
bg.color = "black")+
tm_add_legend(col = c("black", "darkturquoise", "magenta", "purple", "greenyellow","white", "red"),
labels = c("unflagged cells", "cells with 0 total counts", "cells with ctrl/total ratio > 0.1", "outlier for DAPI lower thr.", "outlier for DAPI higher thr.", "outlier for um area lower thr.","outlier for um area higher thr."))
## Warning: The projection of the shape object polygons4 is not known, while it
## seems to be projected.
## Warning: Current projection of shape polygons4 unknown and cannot be
## determined.
## Some legend labels were too wide. These labels have been resized to 0.62, 0.65, 0.63. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
polygons5 <- subset(polygons_sub, subset = fov == 5)
tm_shape(polygons5) +
tm_fill(col= "outlier_fence_color") +
tm_borders(lwd = 0.2, col = "cyan") +
tm_layout(legend.outside = TRUE,
main.title.position = c("center", "top"),
main.title = "",
main.title.fontface = 2,
main.title.size = 1,
inner.margins = c(0, 0, 0, 0),
outer.margins = c(0, 0, 0, 0),
bg.color = "black")+
tm_add_legend(col = c("black", "darkturquoise", "magenta", "purple", "greenyellow","white", "red"),
labels = c("unflagged cells", "cells with 0 total counts", "cells with ctrl/total ratio > 0.1", "outlier for DAPI lower thr.", "outlier for DAPI higher thr.", "outlier for um area lower thr.","outlier for um area higher thr."))
## Warning: The projection of the shape object polygons5 is not known, while it
## seems to be projected.
## Warning: Current projection of shape polygons5 unknown and cannot be
## determined.
## Some legend labels were too wide. These labels have been resized to 0.62, 0.65, 0.63. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
polygons6 <- subset(polygons_sub, subset = fov == 6)
tm_shape(polygons6) +
tm_fill(col= "outlier_fence_color") +
tm_borders(lwd = 0.2, col = "cyan") +
tm_layout(legend.outside = TRUE,
main.title.position = c("center", "top"),
main.title = "",
main.title.fontface = 2,
main.title.size = 1,
inner.margins = c(0, 0, 0, 0),
outer.margins = c(0, 0, 0, 0),
bg.color = "black")+
tm_add_legend(col = c("black", "darkturquoise", "magenta", "purple", "greenyellow","white", "red"),
labels = c("unflagged cells", "cells with 0 total counts", "cells with ctrl/total ratio > 0.1", "outlier for DAPI lower thr.", "outlier for DAPI higher thr.", "outlier for um area lower thr.","outlier for um area higher thr."))
## Warning: The projection of the shape object polygons6 is not known, while it
## seems to be projected.
## Warning: Current projection of shape polygons6 unknown and cannot be
## determined.
## Some legend labels were too wide. These labels have been resized to 0.62, 0.65, 0.63. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
polygons7 <- subset(polygons_sub, subset = fov == 7)
tm_shape(polygons7) +
tm_fill(col= "outlier_fence_color") +
tm_borders(lwd = 0.2, col = "cyan") +
tm_layout(legend.outside = TRUE,
main.title.position = c("center", "top"),
main.title = "",
main.title.fontface = 2,
main.title.size = 1,
inner.margins = c(0, 0, 0, 0),
outer.margins = c(0, 0, 0, 0),
bg.color = "black")+
tm_add_legend(col = c("black", "darkturquoise", "magenta", "purple", "greenyellow","white", "red"),
labels = c("unflagged cells", "cells with 0 total counts", "cells with ctrl/total ratio > 0.1", "outlier for DAPI lower thr.", "outlier for DAPI higher thr.", "outlier for um area lower thr.","outlier for um area higher thr."))
## Warning: The projection of the shape object polygons7 is not known, while it
## seems to be projected.
## Warning: Current projection of shape polygons7 unknown and cannot be
## determined.
## Some legend labels were too wide. These labels have been resized to 0.62, 0.65, 0.63. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
polygons8 <- subset(polygons_sub, subset = fov == 8)
tm_shape(polygons8) +
tm_fill(col= "outlier_fence_color") +
tm_borders(lwd = 0.2, col = "cyan") +
tm_layout(legend.outside = TRUE,
main.title.position = c("center", "top"),
main.title = "",
main.title.fontface = 2,
main.title.size = 1,
inner.margins = c(0, 0, 0, 0),
outer.margins = c(0, 0, 0, 0),
bg.color = "black")+
tm_add_legend(col = c("black", "darkturquoise", "magenta", "purple", "greenyellow","white", "red"),
labels = c("unflagged cells", "cells with 0 total counts", "cells with ctrl/total ratio > 0.1", "outlier for DAPI lower thr.", "outlier for DAPI higher thr.", "outlier for um area lower thr.","outlier for um area higher thr."))
## Warning: The projection of the shape object polygons8 is not known, while it
## seems to be projected.
## Warning: Current projection of shape polygons8 unknown and cannot be
## determined.
## Some legend labels were too wide. These labels have been resized to 0.62, 0.65, 0.63. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
spe$polygons$fscore_outlier <- ifelse(spe$polygons$fscore_outlier == "red", "red", "black")
polygons21 <- subset(spe$polygons, subset = fov == 21)
tm_shape(polygons21) +
tm_fill(col= "fscore_outlier") +
tm_borders(lwd = 0.2, col = "cyan") +
tm_layout(legend.outside = TRUE,
main.title.position = c("center", "top"),
main.title = "",
main.title.fontface = 2,
main.title.size = 1,
inner.margins = c(0, 0, 0, 0),
outer.margins = c(0, 0, 0, 0),
bg.color = "black")
## Warning: The projection of the shape object polygons21 is not known, while it
## seems to be projected.
## Warning: Current projection of shape polygons21 unknown and cannot be
## determined.
polygons21$cell_color <- case_when(polygons21$fov_id == "21_704" ~ "darkturquoise",
polygons21$fov_id == "21_2026" ~ "greenyellow",
polygons21$fov_id == "21_453" ~ "magenta",
polygons21$fscore_outlier == "red" ~ "red",
polygons21$fov_id == "21_158" ~ "white",
TRUE ~ "black")
tm_shape(polygons21) +
tm_fill(col= "cell_color") +
tm_borders(lwd = 0.2, col = "cyan") +
tm_layout(legend.outside = TRUE,
main.title.position = c("center", "top"),
main.title = "Flagged cells",
main.title.size = 0.5,
inner.margins = c(0, 0, 0, 0),
outer.margins = c(0, 0, 0, 0),
bg.color = "black")
## Warning: The projection of the shape object polygons21 is not known, while it
## seems to be projected.
## Warning: Current projection of shape polygons21 unknown and cannot be
## determined.
# cell flagged only for target counts on area and flag score: 21_704, darkturquoise
polygons21[polygons21$fov_id == "21_704",]$flag_score # 0.57
## [1] 0.5730289
polygons21[polygons21$fov_id == "21_704",]$target_area_ratio # 1.07
## [1] 1.075963
abs(polygons21[polygons21$fov_id == "21_704",]$log2AspectRatio) # 0.02
## [1] 0.02856915
# cell flagged only for aspect ratio and flag score: 21_2026, greenyellow
polygons21[polygons21$fov_id == "21_2026",]$flag_score # 0.57
## [1] 0.5700525
polygons21[polygons21$fov_id == "21_2026",]$target_area_ratio # 4
## [1] 4.001505
abs(polygons21[polygons21$fov_id == "21_2026",]$log2AspectRatio) # 0.91
## [1] 0.9183862
# cell flagged for both target counts/area, aspect ratio high and flag score: 21_453, magenta
polygons21[polygons21$fov_id == "21_453",]$flag_score # 0.46
## [1] 0.4675118
polygons21[polygons21$fov_id == "21_453",]$target_area_ratio # 2.35
## [1] 2.35455
abs(polygons21[polygons21$fov_id == "21_453",]$log2AspectRatio) # 0.83
## [1] 0.8365013
# unflagged cell, 21_158, white
polygons21[polygons21$fov_id == "21_158",]$flag_score # 0.96
## [1] 0.968864
polygons21[polygons21$fov_id == "21_158",]$target_area_ratio # 11.6
## [1] 11.60135
abs(polygons21[polygons21$fov_id == "21_158",]$log2AspectRatio) # 0.04
## [1] 0.04264434